Variational determination of approximate bright matter- wave soliton solutions in anisotropic traps 
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We consider the ground state of an attractively-interacting atomic Bose-Einstein condensate in a prolate, 
cylindrically symmetric harmonic trap. If a true quasi-one-dimensional limit is realized, then for sufficiently 
weak axial trapping this ground state takes the form of a bright soliton solution of the nonlinear Schrodinger 
equation. Using analytic variational and highly accurate numerical solutions of the Gross-Pitaevskii equation 
we systematically and quantitatively assess how soliton-like this ground state is, over a wide range of trap 
and interaction strengths. Our analysis reveals that the regime in which the ground state is highly soliton-like 
is significantly restricted, and occurs only for experimentally challenging trap anisotropics. This result, and 
our broader identification of regimes in which the ground state is well-approximated by our simple analytic 
variational solution, are relevant to a range of potential experiments involving attractively-interacting Bose- 
Einstein condensates. 
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I. INTRODUCTION 

Bright soHtons are self-focusing, non-dispersive, particle- 
like solitary waves occurring in integrable systems fll, 0]. 
They behave in a particle-like manner, emerging from mu- 
tual collisions intact except for shifts in their position and rel- 
ative phase. Bright soliton solutions of the one-dimensional 
nonlinear Schrodinger equation (NLSE) can be described an- 
alytically using the inverse scattering technique ^ 01 and 
are well-known in the context of focusing nonlinearities in 
optical fibers fllH]. Bright solitary matter-waves in an at- 
tractively interacting atomic Bose-Einstein condensate (BEC) 
represent an intriguing alternative physical realization |@-[^. 
In a mean-field description an atomic BEC obeys the Gross- 
Pitaevskii equation (GPE) [9], a three-dimensional NLSE. 
While in general non-integrable, in a homogeneous, quasi- 
one-dimensional (quasi- ID) limit the GPE reduces to the one- 
dimensional NLSE, thus supporting bright solitons 1 10-15]. 

Outside the quasi- ID limit the GPE continues to support 
bright solitary matter-waves. These exhibit many soliton-like 
characteristics and have been the subject of much experimen- 
tal lH-d] and theoretical fH-lH] investigation. Both bright 
solitons and bright solitary waves are excellent candidates for 
use in atom interferometry [30], as their coherence, spatial lo- 
caUzation and soliton-like dynamics offer a metrological ad- 
vantage in, e.g., the study of atom-surface interactions ll7l l20ll . 
Towards this end, proposals to phase-coherently split bright 
soUtons and bright solitary waves using a scattering potential 
||27|-|29^ and an internal state interference protocol 1 18], and to 
form soliton molecules ll26ll have been explored in the litera- 
ture. However, while the dynamics and collisions of bright 
solitary waves have been explored in detail and have been 
shown to be soliton-like in three-dimensional (3D) parame- 
ter regimes [16-'T8^, less attention has been directed at the 
question of exactly how soliton-like the ground state of the 
system is. In particular, the experimental feasibility of reach- 
ing the quasi- ID limit of an attractively-interacting BEC, and 
hence obtaining a highly soliton-like ground state, remains 
an area lacking a thorough quantitative exploration. Obtain- 
ing such a ground state, in addition to being interesting in 



its own right, would be highly advantageous in experiments 
seeking to probe quantum effects beyond the mean-field de- 
scription It27r-r29.] , and possibly to exploit the effects of macro- 
scopic quantum superposition to enhance metrological pre- 
cision lUllll. Similar concerns regarding adverse residual 
3D effects in interferometric protocols prompted a recent per- 
turbative study of residual 3D effects in highly anisotropic, 
repulsively-interacting BECs ll33ll . 

The potential instability to collapse of attractively- 
interacting BECs 134-44] is the key obstacle to realizing 
soliton-like behavior in a BEC. Previous studies of bright soli- 
tary wave dynamics, using variational and numerical solutions 
of partially-quasi-lD GPEs fH, El HH H [reductions of 
the GPE to a ID equation which retain some 3D character, 
in contrast to the full quasi- ID limit] and the 3D GPE [16^ 
[Tsi [34I1 . have shown the collapse instability to be associated 
with non-soliton-like behavior. However, previous studies of 
bright solitary wave ground states have focused on identifying 
the critical parameters at which collapse occurs. Approaches 
used in these studies include partially-quasi-lD methods lfl2ll . 
variational methods ll46ll using Gaussian |fToi[34ll47|] and soli- 
ton (sech) [ 34I l48ll ansatzes, perturbative methods [49], and 
numerical solutions to the 3D GPE [3A 15, 43,44,48]. In the 
latter case, the collapse threshold parameters have been exten- 
sively mapped out for a range of trap geometries ll43ll44ll . 

In this paper we use analytic variational and highly accurate 
numerical solutions of the stationary GPE to systematically 
and quantitatively assess how soliton-like the ground state of 
an attractively-interacting BEC in a prolate, cylindrically sym- 
metric harmonic trap is, over a wide regime of trap and inter- 
action strengths. Beginning with previously-considered vari- 
ational ansatzes based on Gaussian ifiol [34l l47ll and soliton 
HUM) profiles, we obtain new, analytic variational solutions 
for the GPE ground state. Comparing the soliton-ansatz vari- 
ational solution to highly accurate numerical solutions of the 
stationary GPE, which we calculate over an extensive param- 
eter space, gives a quantitative measure of how soliton-like 
the ground state is. In the regime where the axial and radial 
trap strengths dominate over the interactions, we show that 
the Gaussian ansatz variational solution gives an excellent ap- 
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proximation to the true ground state for all anisotropies; in 
this regime the ground state is not soliton-like. In the regime 
in which the interactions dominate over the axial, but not the 
radial, trap strength we demonstrate that the soliton-ansatz 
variational solution does approximate the true, highly soliton- 
like ground state. However, we show that the goodness of 
the approximation and the extent of this regime, where it ex- 
ists at all, is highly restricted by the collapse instability; even 
at large anisotropies it occupies a narrow window adjacent to 
the regime where interactions begin to dominate over all trap 
strengths, leading to non-quasi-lD, non-soliton-like solutions 
and, ultimately, collapse. 

Our results have substantial practical value for experiments 
using attractively-interacting BECs; primarily they define the 
challenging experimental regime required to realize a highly 
soliton-like ground state, which would be extremely useful 
to observe quantum effects beyond the mean-field descrip- 
tion such as macroscopic superposition of solitons Il27l - l29|] . 
We note that bright solitary wave experiments to date have 
not reached this regime Secondarily, our quantita- 

tive analysis of a wide parameter space provides a picture 
of the ground state in a wide range of possible attractively- 
interacting BEC experiments. In particular, it indicates the 
regimes in which a full numerical solution of the 3D GPE 
is well-approximated by one of our analytic variational solu- 
tions, which are significantly easier and less time-consuming 
to determine. 

The remainder of the paper is structured as follows: Af- 
ter introducing the most general classical field Hamiltonian 
and stationary GPE in Section [III we begin by discussing the 
quasi- ID limit in Section|lII] In Section lTlI Al we define the di- 
mensionless trap frequency y; in the quasi- ID limit this is the 
only free parameter, and all our results are expressed in terms 
of this quantity. Similarly, our variational ansatzes are moti- 
vated by the limiting behaviors of the solution in the quasi- ID 
case; in this case we define them as Gaussian and soliton pro- 
files, parametrized by their axial lengths. In Sections IIIIBI 
and unci we find, analytically, the energy-minimizing axial 
lengths for each ansatz as a function of y. Comparison of the 
resulting ansatz solutions to highly accurate numerical solu- 
tions of the stationary quasi- ID GPE allows us to determine, 
in the quasi- ID limit, the regimes of low y in which highly 
soliton-like ground states can be realized (Section [III Db . We 
then consider the 3D GPE in Section |IV] The system then has 
a second free parameter in addition to y; we choose this to 
be K, the (dimensionless) trap anisotropy, which is defined in 
Section HVAl In Sections HVEl to lIV E| we define 3D Gaussian 
and soliton ansatzes, adapted from their quasi- ID analogs and 
each parametrized by an axial and a radial length, and find 
the energy-minimizing lengths for each ansatz. In general this 
requires only a very simple numerical procedure, and in the 
limit of a waveguide-like trap can be expressed analytically 
(Section [IV Fl ). In Section HV Gl we compare the ansatz solu- 
tions to highly accurate numerical solutions of the stationary 
3D GPE and, in Section lTV HI assess the potential for realizing 
truly soliton-like ground states. Finally, Section \V\ comprises 
the conclusions. 



II. SYSTEM OVERVIEW 

We consider a BEC of atoms of mass m and (attractive) 
i-wave scattering length a,, < 0, held within a cylindrically 
symmetric, prolate (the radial frequency cj^ is greater than the 
axial frequency w,) harmonic trap. The ground state is de- 
scribed by the stationary Gross-Pitaevskii equation 
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<A(r) = 0, (1) 



where the trapping potential V{r) - m[aj^^x^/2+aj^(y^+z^)/2], 
/I is a real eigenvalue, and the Gross-Pitaevskii wavefunction 
i^(r) is normalized to one. This equation is generated by the 
classical field Hamiltonian (through the functional derivative 
6H[tfr]/6iff* = Aiff) 



dr 



— |V^(r)|2 + V(r)|iA(r)|2 
zm 
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(2) 



This functional of the classical field i/f describes the total en- 
ergy per particle, and the ground state solution minimizes the 
value of this functional. 

When dealing with variational ansatzes for the ground state 
solution, we proceed by analytically minimizing an energy 
functional in the same form as Eq. (|2]i for a given ansatz. 
In contrast, highly accurate numerical ground states are more 
conveniently obtained by solving a stationary GPE of the same 
form as Eq. ([Til. 



III. QUASI-ID LIMIT 
A. Reduction to ID and reseating 

For sufficiently tight radial confinement {a>r » w.t), such 
that the atom-atom interactions are nonetheless essentially 3D 
[as <K (fi/mcL>, y^^] it is conventional 1 10- 15.1 to assume a re- 
duction to a quasi- ID stationary GPE 



2m dx^ 
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ifr(x) = 0. (3) 



Typically i//{r) is taken to be factorized into ifrix) and the radial 
harmonic ground state {maji /nhy^^ exp{-ma)r[y^ + z^]/2fi), 
such that = 2liciJi\as\. Alternative factorizations are also 
possible, which lead to an effective ID equation retaining 
more 3D character than Eq. (O lil2ilT3l[2il[45ll ; similar factor- 
izations have also been introduced for axially rotating BECs 
llsoll and for quasi-2D BECs in oblate traps [51]. In the ab- 
sence of the axial harmonic confining potential {a>^ — > 0), 
there exist exact bright soliton solutions to Eq. (O of the gen- 
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(4) 

where b, = tP' Imgi^N is a length scale characterizing the soli- 
ton's spatial extent, v is the soliton velocity, C is an arbitrary 
displacement, and D is an arbitrary phase. 

This effective ID Gross-Pitaevskii equation contains two 
key length scales: the axial harmonic length a, = (Tilmux)^!^, 
and the soliton length bx- A mathematically convenient way 
to express the single free parameter of Eq. Q is as the square 
of the ratio of these two length scales; 



y = 



(5) 



This parametrization is achieved by working in "soliton 
units"; lengths are expressed in units of b,^ and energies are 
expressed in units of mg^^N' /fP'. This system can be codified 
as ^ = m = ^idA^ - 1, and yields the dimensionless quasi-lD 
GPE 



1 



(6) 



in which y can be interpreted as a dimensionless trap fre- 
quency L15J . The corresponding classical field Hamiltonian 

is 



i^lD[<A] 



dx 



ax 
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(7) 

The choice of y for the single free parameter in the ID GPE 
[Eq. (|6]l] and the classical field Hamiltonian [Eq. (I?)] can be 
most directly pictured as choosing to hold interaction strength 
constant while varying the axial trap strength, parametrized 
by y. Experimentally, however, any of oj,^, w,., a^, and N 
may be varied in order to vary y. In the case y = the ex- 
act ground state solution is a single, stationary bright soliton: 
= sech(x/2)/2. In the following subsections we develop 
analytic variational solutions )p(x) for general j. Compar- 
ing these solutions to highly accurate numerical solutions of 
the quasi- ID GPE then gives a picture of the behavior of the 
ground state with y. Furthermore, these quasi- ID variational 
solutions motivate the later 3D variational solutions and yield 
several mathematical expressions which reappear in the more 
complex 3D calculations. 



B. Variational solution: Gaussian ansatz 



where the variational parameter, Iq, quantifies the axial 
length. In the trap-dominated limit (y — » oo), the true solu- 
tion tends to a Gaussian with Iq-\. Substituting Eq. dSJ into 
Eq. dTji yields (using identities from Appendix [All 
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where Hiu is now expressed as a function of the axial length 
fc- Setting dHiu/dtc - reveals that the variational energy 
described by Eq. (|9]l is minimized when is a positive, real 
solution to the quartic equation 
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The positive, real solution to this quartic is (see solution in 
Appendix IbIi 
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where we have, for notational convenience, defined ;t' to have 
y-dependence such that 
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C. Variational solution: soliton ansatz 



(12) 



Secondly, we consider a sohton ansatz 



(13) 



where the variational parameter, ^s, again quantifies the axial 
length. In the axially un-trapped limit (y — » 0), the true so- 
lution tends to a classical bright soliton, as described by the 
above ansatz with = 1 • The variational energy per particle 
is given by (using identities from AppendixlAli 
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(14) 



which is minimized when 



We first consider the Gaussian variational ansatz 
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Again, this quartic can be solved analytically (see solution in 
Appendix iBt to give the positive, real minimizing value of {$; 



' Equation describes solutions of unit norm. More general soliton solu- 
tions (B/2hl'^)sech{Blx-vt+C]/2h,)e'''^-'-«^^^^ 

have norm B (and eifective mass = B/4), as arise when considering sev- 
eral solitons simultaneously. 
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FIG. 1. Comparison of quasi-lD variational and numerical solutions: 
(a) Energy-minimizing axial lengths {q (Gaussian ansatz, squares) 
and (soliton ansatz, circles) for the quasi-lD GPE. (b) Minimum 
variational energy compared with the numerically calculated ground 
state energy E\\y (black line) for each ansatz: for low y we show 
//iD (solid symbols), which tends to -1/24 as y ^ 0; for high y 
we show H\q = HiD/y (hollow symbols), which tends to 1/2 as 
y ^ oo (//Jq is equal to the energy expressed in the "harmonic 
units," h = m = a>^. = 1). (c) Relative error in the variational energy, 
A = (//id - /?id)//?id- (d) Normalized maximum deformation of the 
best-fitting ansatz wavefunction i/' Ansatz with respect to the numeri- 
cal ground state if/o, Atf/ = max(|i/'Ansatz - i/'ol)/max(i/'o), expressed as 
a percentage. For clarity in (a,b) [(c)], every 16th [20th] datum is 
marked by a symbol. 



ground state i//o(x), and the corresponding ground state en- 
ergy Eiu, uses a pseudospectral method in a basis of sym- 
metric Gauss-Hermite functions; this is a simpHfied version 
of the pseudospectral method used for 3D calculations, which 
is explained in more detail in the next section. Several quan- 
tities are compared in Fig. [Hb-d): the variational minimum 
energies Hio for each ansatz and the numerical ground state 
energy Ei£, are shown in Fig.[Tlb); the relative error between 
//id and Ei£,, defined as A = (//id - £id)/|£idL is shown for 
each ansatz in Fig.lTJc); and the maximum difference between 
the most appropriate ansatz wavefunction (that with lowest A) 
and the numerical ground state wavefunction, expressed as a 
percentage of the maximum value of the numerically exact 
ground state, Ai/r = max(|i/r Ansatz - i/'ol)/max(i/ro) [Fig.fljd)]. 
All the shown computed quantities are insensitive to a dou- 
bling of the numerical basis size from 500 to 1000 states. 

Both the Gaussian and soliton ansatzes provide an excel- 
lent approximation to the exact solutions over a large range 
of y. In the regimes where the relative error in the energy A 
becomes significantly lower than 10 in particular, the dif- 
ference between the ansatz solutions and numerical solutions 
becomes generally indistinguishable from numerical round- 
off error. For the Gaussian ansatz the convergence to this 
regime is noticeably slower than for the soliton ansatz [Fig. 
[lie)]. This effect is a consequence of the parametrization in 
terms of y and the corresponding "soliton units": increasing 
y leads not only to to higher trap strength, but also to higher 
peak densities \il/(x)[^, and hence a stronger nonlinear effect. 

For later comparison to the 3D case, it is useful to define a 
benchmark value of the relative error A that indicates excellent 
agreement between the ansatz and the numerically exact solu- 
tion. Such a definition, however, will vary according to pur- 
pose. As our objectives in this paper relate significantly to the 
shape of the ground state, this forms the basis of our bench- 
mark; a maximum deformation of the wavefunction below 
0.1% of the peak value [as measured by Aip in Fig.fTJd)] cor- 
responds very closely to A < 10"^. Because the relative eiTor 
A saturates to a background value of a; 10"' in regimes where 
the chosen ansatz is inapplicable, a value of A four orders of 
magnitude below this background value thus corresponds to 
an excellent match in shape between the ansatz and the nu- 
merically exact solution. With respect to this benchmark, the 
Gaussian ansatz represents an excellent fit for logio(y) > 1.15, 
while the ground state is highly soliton-like (the soliton ansatz 
represents an excellent fit) for log[Q(7) < -0.95. 



D. Analysis and comparison to ID numerical solutions 

The energy-minimizing axial lengths £g and ^s, defined by 
Eq. ( fTTT i and Eq. ( fTSI ) respectively, are shown as a function of 
y in Fig. \lla). There is no collapse instability in the quasi- 
lD GPE, and solutions are obtained for all (positive, real) y. 
As intended by the chosen forms of the ansatzes, the limiting 
cases are /'g ^ 1 as y ^ oo and ^ 1 as y ^ 0. To eval- 
uate the accuracy of the ansatzes for general y, we compare 
each ansatz with the numerically determined ground state of 
the quasi- ID GPE. The computation of a numerically exact 



IV. BRIGHT SOLITARY WAVE GROUND STATES IN 3D 

A. Rescaling to effective ID soliton units 

We now consider the cylindrically symmetric 3D Gross- 
Pitaevskii equation [Eq. ([T]i]. Compared to the quasi- ID effec- 
tive Gross-Pitaevskii equation [Eq. ©J, three-dimensionality 
introduces an additional relevant length scale, the radial har- 
monic length fl, - {fi/ina>ry^^. We incorporate this into the 
dimensionless trap anisotropy k = Wr/w.t, which forms an ad- 
ditional free parameter Expressed in the same "soliton units" 
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as Eq. Eq. ([T]l becomes 



-W + V(r)-—\,^(r)f-A 

2 KJ 



with corresponding energy functional 



m = 0, (17) 



dv 



-ViA(r) ■ ViA*(r) 



+ y(r)|iA(r)P - — |.A(r) 



, (18) 



where y(r) = y^j^i + ^2(^2 + ^2-)]/2. 

In the following subsections we obtain variational solutions 
for general k and y using ansatzes similar to the Gaussian 
and soliton ansatzes employed in the previous section, with 
an additional variable-width Gaussian radial profile. Contrary 
to the case in the quasi- ID limit, a self-consistent energy- 
minimizing solution for both the axial and radial length pa- 
rameters cannot be expressed entirely analytically. However, 
we reduce the numerical work required to the simultaneous 
solution of two equations, and introduce a straightforward it- 
erative technique to achieve this. We also consider the case 
of a waveguide-like trap (w^ = 0) separately, where an en- 
tirely analytic variational solution exists (Section [IV Fb . Sub- 
sequently, in Section lTVGI we again compare the ansatz solu- 
tions to high-accuracy numerics. 



B. Variational solution: Gaussian ansatz 

We first consider an ansatz composed of Gaussian axial and 
radial profiles. We phrase this as 



(19) 



Here, the first variational parameter, Iq, quantifies the axial 
length of the ansatz in analogy to the quasi- ID case. The re- 
ciprocal of the second variational parameter, fe^', quantifies 
the radial length of the ansatz. In the trap-dominated limit 
(j — > oo) both these lengths approach unity ({^g,^g) — > 1)- 
Substitution of this ansatz into Eq. ( fTSl l yields (using identi- 
ties from Appendix IaIi 



Hioi^ckc) - — 



(20) 



Setting the partial derivatives with respect to both £c and kc 
equal to zero, we deduce that ic must solve the quartic equa- 
tion 



(2ny) 



1/2 



-1=0, 



and that must solve 



{Inyyi^Klc - 1 



1/4 



(21) 



(22) 



From Eq. ( |22] | it follows that we must have {q > 
l/{2nyy^^K to obtain a physically reasonable solution, i.e., a 
real, positive value of kc, consistent with our initial ansatz. 
For a given such value of kc, Eq. (|2TI) is solved (see solution 
in AppendixlBll by 



1/2^2/3 



24/3 (;rr) 



3/2 



1/2 



(23) 



with;^f defined as in Eq. (fT2l l. 



C. Analysis of Gaussian ansatz solution 

Contrary to the quasi- ID limit, minimization of the vari- 
ational energy in 3D requires simultaneous solution of two 
equations for the radial length, k^^, and the axial length, {c- 
These equations are, respectively, Eq. (l22l) and [rearranged 
fromEq. (|2B] 



kn - 
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(24) 



These equations dictate that physical solutions must have 



1 



(Iny) 



1/2, 



<^G < 1, 



(25) 



and hence that y > I /2nK^ must be satisfied in order for phys- 
ical solutions to exist. 

Where solutions exist, they must be found numerically. 
However, a very practical method of numerical solution fol- 
lows from the shape of the dc surface defined by Eq. (l23l l. and 
shown in Fig. Ha), which is a decreasing function of kc for 
all (real, positive) y. The method can be considered graphi- 
cally, in terms of locating the intersection(s) of Eq. (l22b and 
Eq. (|24] |. These curves are shown, for various k, in Fig. |2b- 
d), along with the lower bound from inequality ( |25] ). Below 
a /f-dependent threshold value of y the curves fail to intersect, 
indicating instability of the BEC to collapse. At the thresh- 
old value [dotted curves in Fig. EJb-d)] there is exactly one 
intersection, and above the threshold value [other curves in 
Fig.|2b-d)] there are two intersections. In the latter case the 
higher-^'o intersection, which smoothly deforms to the limit- 
ing case {^g,^g) — ^ 1 as y — » 00, represents the physical, 
minimal-energy variational solution. This solution can be lo- 
cated using a simple "staircase" method: substituting a trial 
value kc, satisfying 1 < < kc, into Eq. |23]produces a trial 
value, ic, satisfying < ^g ^ 1, and subsequently substitut- 
ing this trial value into Eq. |22]produces an iterated trial value, 
satisfying kc < k'^ < kc- Thus, beginning with = 1, 
iteration of this process converges the trial values to the true 
kc and fc- 

The physical solutions to equations i2T[ and (|22] | for differ- 
ent anisotropics k are shown on the {c surface, and projected 
into the ic-y plane, in Fig. |2ja). These solutions are also 
shown as black crosses in the fc-kc plane in Figs. |2b-d). 
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FIG. 2. Energy-minimizing variational parameters for the 3D GPE using a Gaussian ansatz: (a) axial length Iq as a function of the radial length 
and the parameter 7 [Eg. 1231 . Lines show the simultaneous solutions of equations l l22t and l l24t for the axial length £q and radial length 
, for different anisotropies k and values of 7. Projections of these solutions on the j-Cq plane are also shown; here the black line indicates 
the quasi-lD result [from figure[TJa)]. (b-d) Illustration of the intersections of equations J22t [lines with vertical asymptote Cq = l/(2nyy^^K 
shown with fine dashes] and I l24t for various k: the higher-^o intersection, which corresponds to a physical solution for the axial length and 
radial length k^' , can be found using a "staircase" method starting from = I. The numerical solutions obtained this way, and shown by 
points in (a), are shown by crosses in (b-d). The lowest values of 7 plotted in (b-d) are the lowest for which a self-consistent Gaussian ansatz 
solution is found. 



where they form a line connecting the physical-solution in- 
tersections of Eq. (I22I 1 and Eq. ( l24l i for the various y shown. 
In Fig. |3a) the collapse instability is manifest as a rapid rise 
in kc — corresponding to a decrease in radial extent — and 
fall in £g — corresponding to a decrease in axial extent — 
just above a /c-dependent threshold value of y. There are no 
self-consistent solutions for these quantities below this col- 
lapse threshold. For increasing anisotropies a:, this collapse 
threshold occurs at lower values of y. For the highest two val- 
ues of K considered the collapse threshold lies in the regime 
where (q is already approaching 0; our analysis of the Gaus- 
sian ansatz in the quasi-lD limit indicates that the 3D Gaus- 
sian ansatz will be a poor approximation to the true solution 



in this regime. Importantly, for y above the collapse thresh- 
old the projected curves for each anisotropy agree well with 
the Gaussian ansatz in the quasi- ID GPE, suggesting that the 
Gaussian ansatz gives a good approximation to the true solu- 
tion here. 



D. Variational solution: soliton ansatz 

Secondly, we consider a soliton ansatz composed of a axial 
sech profile and a radial Gaussian profile. We phrase this as 

1/2 1/2^ ^ , 

"^^•■^ " TTTW^"'*^^' ""'^'sech(^/2^s). (26) 
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FIG. 3. Energy-minimizing variational parameters for the 3D GPE using a soliton ansatz: (a) axial length as a function of the radial length 
fcj ' and the parameter y [Eg. 1301 . Lines show the simultaneous solutions of equations J33l > and J34t for the axial length and radial length fcg ' , 
for different anisotropies k and values of y. Projections of these solutions on the y~£s plane are also shown; here the black line indicates the 
quasi-lD result [from figurelUa)]. (b-d) Illustration of the intersections of equations ([33} [lines with vertical asymptote = (n/3)/(2nyy'^K 
shown with fine dashes] and i34l for various k: the higher-fs intersection, which corresponds to a physical solution for the axial length and 
radial length can be found using a "staircase" method starting from kg = I. The numerical solutions obtained this way, and shown by 
points in (a), are shown by crosses in (b-d). The lowest values of y plotted in (b-d) are the lowest for which a self-consistent soliton ansatz 
solution is found. 



As with the 3D Gaussian ansatz, the first variational parame- 
ter, fa, quantifies the axial length of the ansatz and the recip- 
rocal of the second variational parameter, fe^', quantifies its 
radial length. In the quasi- ID limit both lengths consequently 
approach unity ({^G, ^g) — > !)■ Substituting this ansatz into 
Eq. dTSl ) yields (using identities from AppendixlAl) 



2 2 



3Kkl 



3k 



(27) 



Once again, setting partial derivatives with respect to both /'s 
and equal to zero allows us to deduce that: 



1 



S /^„2y2 4„2y2 

and that must solve 



0, 



6Ky£s - 1 



1/4 



(28) 



(29) 



From Eq. ( |29] | it follows that we must have > l/6Ky 
to obtain a physically reasonable solution, i.e., a real, positive 
value of ks, consistent with our initial ansatz. For a given such 
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value of ^s, Eq. (1281 ) is solved (see solution in AppendixlBt by 



2ii/6(;ry)2/3 
with;^' defined as in Eq. ( fT2] ) 



3/2 



- 1 



1/2 



(30) 



E. Analysis of soliton ansatz solution 

As in the case of the Gaussian ansatz, minimization of the 
variational energy in 3D requires the simultaneous solution 
of equations for the radial length and the axial length ^s- 
These equations are, respectively, Eq. (|29l ) and [rearranged 
from Eq. ^] 



1/2 



(31) 



These equations dictate that physical solutions must have 

6Ky ^ (2;ry)'/2 

and hence that y > (ji jj)^ jlnK^ must be satisfied in order for 
physical solutions to exist. These equations and constraints 
can be further simplified by casting them in terms of {'^ - 
(Inyyi^i^; this yields two equations. 



1/4 



and 



(27rr)i/2 



and an inequality. 



n/3 



1/2 



(Iny) 



1/2, 



< 4 < 1, 



(33) 



(34) 



(35) 



which are extremely similar to those encountered in the case 
of the Gaussian ansatz. The numerical solution of these equa- 
tions for the physical solution, which can only exist when 
y > {nl3fl2nK^, follows the same procedure as used for the 
Gaussian ansatz. 

Variational-energy-minimizing solutions to the soliton 
ansatz equations for different anisotropics k are shown in Fig. 
[3j these are shown superimposed on the surface and pro- 
jected into the ^s-y plane in Fig.^a), and alongside equations 
( l28T l and ^ and inequaUty in Fig.Ob-d). The collapse 
instability is even more evident in the soliton ansatz than in the 
Gaussian ansatz, since it occurs in a region with a larger back- 
ground value of ^s- Once again, the collapse is manifest as a 
rapid rise in and drop in — corresponding to both axial 
and radial contraction of the solution — immediately prior to 
a /c-dependent threshold value of y. Below the threshold, no 



self-consistent solutions exist. For increasing anisotropics k, 
this collapse threshold again occurs at lower values of y. In 
contrast to the case of the Gaussian ansatz, however, the col- 
lapse instability precludes solutions in exactly the limit where 
one expects the soliton ansatz to be accurate (y — > 0). This 
property of the collapse instability severely restricts the possi- 
bility of observing highly bright-soliton-like ground states in 
3D. The solution curves in Fig. Oa) illustrate that this effect 
is worst for low trap anisotropics k, but is to some extent mit- 
igated for higher k. However, a full comparison with numeri- 
cally exact solutions is necessary to quantify these effects; we 
undertake such a comparison in Section lTV Gl 



F. Variational solution: waveguide configuration 

In broad experimental terms, the collapse instability sets 
a maximum value for the ratio of interaction strength to trap 
strength (equivalent to a minimum value of y) which increases 
(and hence the minimum value of y decreases) with the trap 
anisotropy k. In the context of atomic BEC experiments 
one would typically think of controlling the interaction-trap 
strength ratio by varying either \as\ or while holding w,- and 
w V constant; in this situation the collapse instability places a 
trap-anisotropy-dependent upper limit on the product |a.5|A^. 
However, the minimum value of y does not increase without 
limit in the trap anisotropy k: In an experiment one can, in 
principle, remove all axial trapping to create a waveguide-like 
configuration; in this case Wv = and the trap anisotropy k — > 
oo, while the parameter y — > 0. In this limit a reparametriza- 
tion is necessary, and only needs to be performed for the soli- 
ton ansatz, which is clearly more appropriate in this context. 

Elimination of the axial trap eliminates one of the two free 
parameters of the 3D GPE [Eq. (fT7]i1. The remaining free 
parameter is Y - yK - {arl2\as\Nf', where - {H/mcL>, y^^ is 
the radial harmonic oscillator length scale. The soliton ansatz 
may be re-written in terms of F as 



sech(x/2/'s). 



(36) 



Substituting this into Eq. ( fTSI ) with W j - yields (using iden- 
tities from AppendixlAb. 



( 1 



^ + — + 



24-il 12 £s ' 2 ' 2kl^ 



(37) 



from which we deduce that the energy-minimizing variational 
parameters satisfy 



1 

k\ 



and 



/_6Rs_ 
Wh - 1 



1/4 



(38) 



(39) 



Contrary to the more general 3D case, an analytic simultane- 
ous solution of Eq. ( |38] ) and Eq. ( [39] ) exists when satisfies 
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c 
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60 

60 
C 




1 2 

logio(r) 

FIG. 4. Comparison of 3D variational and numerical solutions in 
a waveguide configuration (aij = 0: (a) Energy-minimizing ax- 
ial length and radial length fcj' for the soliton ansatz. Solu- 
tions, given by Eq. | |4U . exist for all T = kj > 3''^/4. (b) Rela- 
tive error in the minimum variational energy of the soliton ansatz, 
A = (HjD - Eij))/E2j), where £30 is the numerically determined 
ground state energy. 



the depressed cubic equation 



(40) 



Using the general solution for a depressed cubic equation from 
Appendix IbI one finds that the physical root (with real, posi- 
tive (s satisfying the limit /"s — > 1 as T ^ 00) is given by 



i2r 33/2r\i6 



1 



1/2 



1/3 



i2r 33/2r \ 16 



1/2 



1/3 



(41) 



Consequently, solutions only exist for T > 3'^^/4, as shown in 
Fig. Ha). 



G. Comparison to 3D numerical solutions 

The variational energy-minimizing axial lengths £g and 
are shown as functions of y in Fig. |5ja) for the general 3D 
case; for the waveguide limit both axial and radial lengths 
is and feg' are shown as functions of F in Fig. |5a). As in 
the quasi- ID case, we quantitatively evaluate the accuracy of 
the ansatz solutions for general y (F) by comparing the varia- 
tional minimum energy H^d with the numerically determined 
ground state energy £30. We calculate E^u using a pseu- 
dospectral method in a basis of optimally-scaled harmonic 
oscillator eigenstates; this is formed from a tensor product 
of symmetric Gauss-Hermite functions (axial direction) and 
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FIG. 5. Comparison of 3D variational and numerical solutions: (a) 
Energy-minimizing axial lengths {q (Gaussian ansatz, solid symbols) 
and Cs (soliton ansatz, hollow symbols), (b) Scaled variational ener- 
gies //jp = KHij)/y(K + 1/2) (a similarly scaled ground state en- 
ergy = KEiolyiK + 1/2) tends to 1 in the limit y ^ 00 for 
all anisotropics k) compared with the numerically calculated ground 
state energies £30 (black dots). (c,d) Normalized relative error in the 
variational energy A = (//3d-£'3d)/£'3d for the Gaussian (c) and soli- 
ton (d) ansatzes. For clarity every 4th datum is marked by a symbol 
in (a-d). 



generalized Laguerre functions (radial direction). The ansatz 
with the lowest variational energy is used both to optimize the 
scaling of the basis functions and as an initial estimate for the 
solution. Expanding the stationary 3D GPE in such a basis 
produces a system of nonlinear equations which are solved it- 
eratively using a modified Newton method. A similar method 
was used to solve a similar cylindrically symmetric, stationary 
3D GPE, with repulsive interactions, in Ref. i52ll . 

As in the quasi- ID case, we compare several quantities be- 
tween the ansatz and numerical solutions. Fig.|5|b) shows the 
scaled energy H'^^ - {HT,Y)ly)l{\ + I 12k) in the general 3D 
case. This scaling is such that — which is defined anal- 
ogously to with respect to £30 — tends to 1 as y — » 00. 
Figs. |3c) and (d) show the relative error in the variational 
minimum energy A = {H^d - E3d)/E3d for the Gaussian and 



10 



soliton ansatzes, respectively. The same quantity A is shown 
for the waveguide limit in Fig. |4|b). All quantities shown in 
Figs. |5] and m are computed using between 2000 and 12000 
basis states (/c-dependent) and are insensitive to a doubling of 
the number of basis states. 

In the general 3D case, a close inspection of Fig. lUb-d) 
is necessary to reveal the overall relation between the ansatz 
solutions and the numerically obtained ground state. In the 
high-y limit Fig. |3b) shows that both the Gaussian varia- 
tional energies (solid symbols) and the ground state energy 
£30 (black dots) approach 1 as y — » 00, whereas the soliton 
ansatz energies (hollow symbols) tend to higher energies. This 
corresponds to the actual ground state most closely matching 
the Gaussian ansatz in this limit, as one would expect. In- 
deed, the relative error in variational energy. A, for the Gaus- 
sian ansatz [Fig.|3c)] continues to drop exponentially with j 
for all anisotropies a:, making it possible to find regimes of y 
where the Gaussian ansatz gives an excellent approximation 
to the true ground state. 

In the opposite, low-y limit, collapse occurs at a k- 
dependent value of y; this corresponds to the points in Fig. 
|5ja-d) where solution curves abruptly cease. Prior to col- 
lapse (at higher values of y) the relation between the Gaussian 
ansatz, the soliton ansatz, and the actual ground state is highly 
dependent on the trap anisotropy k [Fig. |5lb)]. In the case 
of a spherically symmetric trap, where the anisotropy k - \, 
the soliton ansatz variational energy is never closer to the true 
ground state energy £30 than the Gaussian ansatz variational 
energy. A regime of soliton-like ground states consequently 
cannot exist at this low anisotropy; as the soliton ansatz is 
intrinsically asymmetric, this is to be expected. For higher 
anisotropies, the soliton ansatz energy is closer to £30 than 
the Gaussian ansatz energy in a small regime prior to collapse. 
Exactly how soliton-like the ground state is in this regime can 
be quantitatively assessed using the relative error A. This is 
shown for the soliton ansatz in [Fig. |3d)]. For each k the 
"background" value of A in the limit -y ^ 00 is different; this 
effect is due to the decreasing size of the axial part of the en- 
ergy with respect to the radial part for increasing y. In the op- 
posite, low-y, limit A increases sharply close to the collapse 
point as the ground state wavefunction rapidly contracts. The 
maximum extent to which A decreases from its high-y limit, 
before this increase due to collapse-related contraction at low 
y, quantifies how soliton-like the ground state becomes in this 
regime. Even for the highest anisotropy shown, k - 256, the 
regime of y over which A drops below its background value 
is rather narrow, and the actual drop in A is only one order 
of magnitude. Compared to benchmark of Section ITlIDI this 
indicates that the true ground state remains considerably de- 
formed with respect to the soliton ansatz. The minimum error 
in the soliton ansatz energy does, however, improve with in- 
creasing anisotropy k. Excellent agreement can be achieved 
in the waveguide limit (k 00): Fig. |4] shows that excellent 
agreement, with respect to the benchmark figure of Section 
imp] can be obtained for F > 10^^^ 



H. Discussion 

A physical interpretation of the above results follows from 
considering two conditions that must be satisfied in order to 
realize a soliton-like ground state; (1) the radial profile should 
be "frozen" to a Gaussian, thus realizing a quasi- ID limit; and 
(2) interactions should dominate over the axial trapping. On 
first inspection these conditions seem mutually compatible, 
and satisfiable simply by increasing the radial trap frequency 
(jjr with other parameters held constant. However, condition 
(1) can only be satisfied if the maximum density remains low 
enough to avoid any deformation of the radial profile due to 
the collapse instability. Increasing w,- leads to exactly such de- 
formation, and ultimately to collapse, as it has the secondary 
effect of strongly increasing the density. This strong increase 
in density with w,. is particular to the case of attractive interac- 
tions. Increasing w,. in a repulsively-interacting BEC likewise 
acts to increase the density, but this increase is counteracted 
by the interactions; these act to reduce the density, and cause 
the BEC to expand axially. In the attractively-interacting case 
the response of the interactions is the opposite: increasing w, 
leads to axial contraction of the BEC. Consequently condi- 
tion (1) is far harder to satisfy for an attractively-interacting 
BEC than a repulsively-interacting one. Responding to this 
problem simply by reducing the interaction strength (either 
through Ifl^l or A^) leads to violation of condition (2). The 
nature of the problem is made particularly clear by consid- 
ering the waveguide limit; here condition (2) is automati- 
cally satisfied {lj^ — 0). This makes it possible to achieve 
a highly soliton-like ground state by satisfying condition (1) 
alone. However, such a ground state is achieved by lowering 

1 /2 

the product and thus by progressing towards the 

limit of extreme diluteness. 

This physical behavior of the system presents considerable 
challenges for experiments aiming to realize a highly soliton- 
like ground state. In essence, the most desirable configuration 
is to have extremely high anisotropies k, while keeping w,- as 
low as possible. Realizing such a configuration through ex- 
tremely low, or zero, axial trap frequencies is problematic: 
such frequencies are hard to set precisely experimentally as 
they require a very smooth potential to be generated, poten- 
tially over a considerable length. Furthermore, in the case 
dLt^ = the mean-field approximation ceases to be valid for an 
attractively-interacting BEC; the true wavefunction should be 
translationally invariant in this case, but the mean-field solu- 
tion breaks this symmetry [53]. Even for very low but non- 
zero the mean-field approximation can lose validity due to 
the extreme diluteness of the BEC, and the energy gap from 
the ground state to states with excited axial modes can be- 
come low enough to cause significant population of the ex- 
cited states at experimentally feasible temperatures. 

It is informative to consider the parameters used in bright 
solitary wave experiments to date [6- 8]. None of these aimed 
to realize highly soliton-like ground states in the sense consid- 
ered here. However, they nonetheless indicate regimes which 
have proved to be experimentally accessible and offer a guide 
to future possibilities. All have operated outside the regime 
of highly soliton-like ground states; direct comparison of the 
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experiments of Refs. and f§] with our results reveals that 
K is too small in these experiments (a- ^ 11 and k ~ 3 re- 
spectively) to achieve a highly soli ton-like ground state. The 
experiment of Ref. featured an expulsive axial potential, 
which does not yield a value of k suitable for direct compar- 
ison with our results. However, it is possible to assume the 
waveguide limit cj^ = in each experiment and compare the 
values of F with our results: in each case T < 1, outside the 
regime of highly soliton-like ground states. Thus, experiments 
with weaker traps and lower densities than previously realized 
with attractive condensates appear to be necessary in order to 
achieve a highly soliton-like ground state. 



esting direction for future work. 
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Appendix A: Useful integrals 



V. CONCLUSIONS 

In this paper we considered attractively-interacting atomic 
BECs in cylindrically symmetric, prolate harmonic traps, and 
introduced variational ansatzes, based on Gaussian and bright- 
soliton profiles, for the GPE ground state. We compared new, 
analytic variational solutions based on these ansatzes with 
highly accurate numerical solutions of the GPE over an ex- 
tensive parameter space, and hence determined how soliton- 
like the ground state is. Initially assuming the quasi-lD limit 
to be valid, we showed that the true solution to the GPE is 
(not) soliton Uke when interactions do (not) dominate over 
the trap strength. In 3D, this picture is complicated by the 
collapse instability; in the regime where all trap strengths 
dominate over the interactions a Gaussian variational ansatz 
gives an excellent approximation to the true, and non-soliton- 
like ground state. In contrast to the quasi- ID limit, however, 
we have shown that the regime in which the ground state is 
truly soliton-like (well approximated by a soliton variational 
ansatz) is either non-existent, or highly restricted, depending 
on the trap anisotropy. For low anisotropics, as one raises the 
strength of the interactions such that they approach and ex- 
ceed the strength of the axial trap the true ground state ceases 
to be well-described by a Gaussian variational ansatz, but does 
not become well-described by a soliton variational ansatz be- 
fore the interaction strength also exceeds the radial trapping 
strength, leading to collapse. Only by raising the anisotropy 
significantly can one open a parameter window in which the 
true ground-state becomes soliton-like before the interaction 
strength is sufficient to cause collapse. 

Our results describe the nature of the ground state over a 
wide parameter regime, and offer a straightforward, accurate 
approximation to the full 3D GPE solution in many cases. 
Our results are particularly relevant for experiments using 
attractively-interacting condensates as they identify the po- 
tentially challenging parameter regime required to observe a 
truly soliton-like ground state, which would be an advanta- 
geous regime for experiments seeking to explore and exploit 
beyond-mean field effects such as a macroscopic superposi- 
tion of bright solitons. Given that previous studies have shown 
that the dynamics and collisions of bright solitary waves can 
be soliton-like over a much wider parameter regime than our 
approach reveals the ground state to be, extending the varia- 
tional approach used here to dynamical situations is an inter- 



Considering a Gaussian ansatz to be proportional to e , 
for completeness we reprise the following sequence of well- 
known integral identities, all of which are necessary to deter- 
mine the corresponding variational energy functional: 



•J —Ck 



dxe 



-2k- A 



/^^ r 

1 d r°° 



dxe 



-2k'. 



^, (Al) 
2k ' 



(A2) 



4k^ 



dxi—e-'''-''] =4k* dxx^e-'''-'' ^k^/^. (A3) 
oo \ox j 

Comparable integral identities exist when considering an 
ansatz proportional to sech(A;jic). Thus: 



f 



dxsech (kx) 



tanh(A:x) 



(A4) 



f 

Xco 
oo 



dxsech (kx) 



— sech(fex) 
ox 



{sech (kx) + 2)tanh(fcx) 



3k 



4 

3k' 



(ixtanh (A;jic)sech (jc) 



k ?k 
=1 [tanh(^x)]-„ = -, 



(A5) 



(A6) 



all of which are necessary to determine the energy of a stan- 
dard bright soliton solution to the nonlinear Schrodinger equa- 
tion. However, we also require a contribution arising from the 
existence of an external harmonic confining potential. Hence, 
we determine 



dxx^sech^(kx) -2 I dxx^sech^(kx) 
oo Jo 



^[Li2 [-e-^''-') + kx[kxtmh(kx) 



-kx-2ln[l+e-^''')]]'^ (A7) 
= ^[Li2(0)-Li2(-l)] 

where Liy(x) = 2,7=1 is a polylogarithm, and 

-Liv(-l) = T](y), the Dirichlet rj function, with 7/(2) = l\2. 
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Appendix B: Solution to the quartic equations 

We require a general solution to a quartic in I of the form 



r + - c = 0, 



(Bl) 



where h and c are positive real constants, and I must also take 
positive real values to be physically meaningful. This can be 
rephrased as the product of two quadratics in I: 



= 0, (B2) 



so long as {b^/a^ - a'^)/4 - c. Hence, a, which remains to be 
determined, must solve + Aca^ — b^. 

Defining ^ = a^, the problem of determining a reduces to 
finding values of ^ to solve the depressed cubic equation 



f + 4c^ -b^ ^Q. 



(B3) 



Defining 



A=\\—+\— + , B 

V 2 V 4 27 ' 



^b^ 64^3 



the three roots of Eq. ( IB3b are given by: 
^1 =A + B, 

^2 = - (A + B)I2 + i V3(A - B)/2, 
^3 = - (A + B)/2 - / V3(A - B)I2. 



(B4) 

(B5) 
(B6) 
(B7) 



Any one of these will solve Eq. (IB3l l. however we choose ^\ ; 
as b and c are assumed positive real, ^\ is also conveniently 
guaranteed positive real. 

Substituting in a = we can apply the quadratic for- 
mula to both the factors (enclosed in square brackets) on the 



left hand side of Eq. 



This reveals the four roots to be 





V ^' 


+ 2bl 
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V ^' 


+2^7/ v?r 
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-2bl^[ri 
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-2blyff, 




2 





(B8) 
(B9) 
(BIO) 
(BID 



Recalling that and are positive real, ^3 and {4 are clearly 
complex, and therefore not of interest to us. Noting that 



^\^A^ + B^ + 3AB(A + B)^b-- 4c^i 



(B12) 



we can see that + B^ = b^ > , hence Ab^ > and thus 
2b/ > ^1 • Roots i 1 and {2 are therefore real, but {2 is 
guaranteed negative. However, fromEq. (IB12I) it also follows 
that 



b>^iyffi^2b/ylfi>2^i^2b/yffi-^i 



(B13) 



Hence ^1 is guaranteed positive real, and is the only solution 
of interest. 

Thus, the single positive real root of Eq. (IBll) is 



1/3 



27/6 



3/2 



1/2 



(B14) 



with 



1 + 



(c/3) 



3 -,1/2^ 1/3 



(b/AY 



+ 1- 



1 + 



(c/3) 



3 -,1/2^ 1/3 



ib/4r 



(B15) 

and where values of all fractional powers are taken to be real, 
and positive when a positive root exists. 
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